knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(raster)
library(sf)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))
Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.
Included species are:
spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE) %>%
select(species) %>%
mutate(fe = TRUE) %>%
distinct()
spp_am <- get_am_spp_info() %>%
filter(occur_cells >= 10) %>%
select(species = sciname) %>%
mutate(am_mapped = TRUE) %>%
distinct()
spp_iucn <- read_csv(here('_data/iucn_spp/iucn_to_worms_match.csv'), show_col_types = FALSE) %>%
rename(species = worms_name, iucn_mapped = mapped)
spp_vuln <- get_spp_vuln() %>%
select(species) %>%
mutate(vuln = TRUE) %>%
distinct()
spp_worms <- assemble_worms() %>%
select(species) %>%
mutate(worms = TRUE) %>%
distinct()
all_spp <- spp_worms %>%
full_join(spp_vuln, by = 'species') %>%
full_join(spp_am, by = 'species') %>%
full_join(spp_iucn, by = 'species') %>%
full_join(spp_fe, by = 'species')
spp_for_fe_metrics <- all_spp %>%
### names in WoRMS
filter(worms) %>%
### only those spp in vulnerability project?
# filter(vuln) %>%
### mapped in AquaMaps and/or IUCN
filter(am_mapped | iucn_mapped) %>%
### functional traits available
filter(fe)
# worms_mismatch <- all_spp %>%
# filter(is.na(worms)) %>%
# filter()
#
# ### all AquaMaps spp
# am_check <- get_am_spp_info() %>%
# filter(sciname %in% worms_mismatch$species)
# table(am_check %>% select(kingdom, phylum))
# table(am_check %>% filter(kingdom == 'animalia') %>% select(class, kingdom))
Here we use the IUCN species range maps rasterized to 10 km Mollweide, and AquaMaps species distributions reprojected to 10 km Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥ 10.
Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species (separately for IUCN species and AquaMaps species) for later calculation into functional richness, redundancy, vulnerability.
iucn_fe_summary_file <- here_anx('func_entities', 'fe_cell_summary_iucn.csv')
# unlink(iucn_fe_summary_file)
if(!file.exists(iucn_fe_summary_file)) {
spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE)
iucn_to_map <- spp_for_fe_metrics %>%
left_join(spp_fe, by = 'species') %>%
filter(iucn_mapped) %>%
select(fe_id, iucn_sid) %>%
distinct() %>%
arrange(iucn_sid)
iucn_map_fstem <- here_anx('spp_maps_mol', 'iucn_spp_mol_%s.csv')
iucn_map_fs <- sprintf(iucn_map_fstem, iucn_to_map$iucn_sid)
message('Assembling IUCN species maps and binding FEs...')
iucn_maps <- parallel::mclapply(iucn_map_fs, mc.cores = 40,
FUN = function(f) {
data.table::fread(f) %>%
filter(presence != 5) %>%
select(-presence)
}) %>%
setNames(iucn_to_map$fe_id) %>%
data.table::rbindlist(idcol = 'fe_id') %>%
mutate(fe_id = as.integer(fe_id))
iucn_fes <- iucn_to_map$fe_id %>% unique() %>% sort()
### summarize was having a hard time with this massive df, so
### try breaking it down by FE and parallel processing it
message('Summarizing IUCN Functional Richness map...')
iucn_fe_summary <- parallel::mclapply(iucn_fes, mc.cores = 4,
FUN = function(fe) { ### fe <- iucn_fes[1]
df <- iucn_maps %>%
filter(fe_id == fe)
message('Processing fe = ', fe, '... ', nrow(df), ' species-cell combos')
df_out <- df %>%
group_by(cell_id) %>%
summarize(n_spp_fe = n(),
fe_id = fe,
.groups = 'drop')
return(df_out)
}) %>%
data.table::rbindlist() %>%
mutate(src = 'iucn')
rm('iucn_maps')
write_csv(iucn_fe_summary, iucn_fe_summary_file)
}
am_fe_summary_file <- here_anx('func_entities', 'fe_cell_summary_am.csv')
# unlink(am_fe_summary_file)
if(!file.exists(am_fe_summary_file)) {
spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE)
am_to_map <- spp_for_fe_metrics %>%
left_join(spp_fe, by = 'species') %>%
filter(am_mapped) %>%
filter(!iucn_mapped | is.na(iucn_mapped)) %>%
select(fe_id, species) %>%
distinct() %>%
arrange(species)
am_map_fstem <- here_anx('spp_maps_mol', 'am_spp_mol_%s.csv')
am_map_fs <- sprintf(am_map_fstem, str_replace_all(am_to_map$species, ' ', '_'))
message('Assembling Aquamaps species maps and binding FEs...')
am_maps <- parallel::mclapply(am_map_fs, mc.cores = 40,
FUN = function(f) {
data.table::fread(f) %>%
filter(prob >= 0.5) %>%
select(-prob)
}) %>%
setNames(am_to_map$fe_id) %>%
data.table::rbindlist(idcol = 'fe_id') %>%
mutate(fe_id = as.integer(fe_id))
am_fes <- am_to_map$fe_id %>% unique() %>% sort()
### summarize was having a hard time with this massive df, so
### try breaking it down by FE and parallel processing it
message('Summarizing Aquamaps Functional Richness map...')
# am_fe_summary <- parallel::mclapply(am_fes, mc.cores = 2,
am_fe_summary <- lapply(am_fes, ### for memory issues, just use a single core... ugh
FUN = function(fe) { ### fe <- am_fes[1]
df <- am_maps %>%
filter(fe_id == fe)
message('Processing fe = ', fe, '... ', nrow(df), ' species-cell combos')
df_out <- df %>%
group_by(cell_id) %>%
summarize(n_spp_fe = n(),
fe_id = fe,
.groups = 'drop')
return(df_out)
}) %>%
data.table::rbindlist() %>%
mutate(src = 'am')
rm('am_maps')
write_csv(am_fe_summary, am_fe_summary_file)
}
tot_fe_summary_file <- here_anx('func_entities/fe_cell_summary_total.csv')
### unlink(tot_fe_summary_file)
if(!file.exists(tot_fe_summary_file)) {
am_fe_summary <- data.table::fread(am_fe_summary_file)
iucn_fe_summary <- data.table::fread(iucn_fe_summary_file)
fe_sum_maps <- data.table::rbindlist(list(am_fe_summary, iucn_fe_summary))
tot_fes <- fe_sum_maps$fe_id %>% unique() %>% sort()
tot_fe_summary <- parallel::mclapply(tot_fes, mc.cores = 12,
FUN = function(fe) { ### fe <- tot_fes[1]
message('Summarizing across AM and IUCN for FE #', fe)
df <- fe_sum_maps %>%
filter(fe_id == fe) %>%
group_by(cell_id) %>%
summarize(n_spp_fe = sum(n_spp_fe, na.rm = TRUE),
fe_id = fe,
.groups = 'drop')
return(df)
}) %>%
data.table::rbindlist()
write_csv(tot_fe_summary, tot_fe_summary_file)
}
Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.
This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')
### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
# system.time({
fe_sum_total <- data.table::fread(tot_fe_summary_file)
### there are 6.5 million cells. Chop into 500k instances and mclapply to summarize
# cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
chunk_size <- 500000
n_chunks <- ceiling(ncell(raster::raster(here('_spatial/ocean_area_mol.tif'))) / chunk_size)
# system.time({
message('Calculating species richness map...')
nspp_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
nspp_sum <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
group_by(cell_id) %>%
summarize(n_spp = sum(n_spp_fe, na.rm = TRUE),
n_fe = n_distinct(fe_id),
.groups = 'drop')
}) %>%
data.table::rbindlist()
# }) ### 30 sec
# system.time({
message('Calculating functional vulnerability map...')
f_vuln_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_vuln <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
group_by(cell_id) %>%
summarize(n_fe = first(n_fe), n_spp = first(n_spp),
f_vuln = sum(n_spp_fe == 1) / n_fe,
f_wvuln = mean(0.5^(n_spp_fe - 1)),
f_red = mean(n_spp_fe),
.groups = 'drop')
}) %>%
data.table::rbindlist()
# }) ### 30 sec
# system.time({
message('Calculating functional overredundancy map...')
f_overred_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_overred <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
mutate(f_over = ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)) %>%
group_by(cell_id) %>%
summarize(f_overred = sum(f_over) / sum(n_spp_fe),
.groups = 'drop')
}) %>%
data.table::rbindlist()
# }) ### 7 sec
message('Combining functional diversity metrics maps...')
fe_metrics_df <- f_vuln_df %>%
oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full') %>%
mutate(across(where(is.numeric), .fns = ~round(.x, 6)))
### rounding to 6 dec places saves a lot of storage for tifs
# })
# ggplot(fe_metrics_df, aes(x = f_vuln, y = f_wvuln)) +
# geom_point(alpha = .3) +
# geom_abline(slope = 1, intercept = 0)
write_csv(fe_metrics_df, fe_metrics_file)
}
Create and save rasters of functional diversity metrics
fe_metrics_df <- data.table::fread(fe_metrics_file)
n_spp_rast <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')
writeRaster(n_spp_rast, here('_output/func_entities/n_spp.tif'),
overwrite = TRUE)
writeRaster(n_fe_rast, here('_output/func_entities/n_fe.tif'),
overwrite = TRUE)
writeRaster(f_vuln_rast, here('_output/func_entities/f_vuln.tif'),
overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'),
overwrite = TRUE)
writeRaster(f_red_rast, here('_output/func_entities/f_red.tif'),
overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'),
overwrite = TRUE)
Examine distributions of each raster.
Plot each raster.